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Abstract 



> 

We study traffic congestion by analyzing a one dimensional traffic flow model. 

o" 

\ Developing an asymptotic method to investigate the long time behavior near 

| a critical point, we derive the modified Korteweg-de Vries (mKdV) equation 

as the lowest order model. There are an infinite number of kink solitons 
to the mKdV equation, while it has been found by numerical simulations 
that the kink pattern arising in traffic congestion is uniquely determined ir- 
respective of initial conditions. In order to resolve this selection problem, we 
consider higher order corrections to the mKdV equation and find that there 
is a kink soliton which can deform continuously with the perturbation rep- 
resented by the addition of these corrections. With numerical confirmation, 
we show that this continuously deformable kink soliton characterizes traffic 
congestion. We also discuss the relationship between traffic congestion and 

slugging phenomenon in granular flow. 
46.10.+Z 47.54. +r 47.20.Ky 
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I. INTRODUCTION 



Cooperative behavior in dissipative systems composed of discrete elements has attracted 
considerable interest. In particular, granular flow exhibited by powder assemblies under 
non-equilibrium conditions |T]-f|], traffic flow and collective motion of living organisms 

such as birds, fish, insects, bacteria and motional protein have been investigated extensively 
PHT2"|, and their spatio-temporal patterns have been reproduced successfully by computer 
simulations of both discrete and continuous models [|T3| -|2Ti|| . These successes suggest that at 
least some of these phenomena can be categorized into one universality class. 

As one example, we find strong resemblance between slugging phenomenon in granular 
flow and congestion in traffic flow. The problem we address here is to uncover the com- 
mon structure behind such apparently similar phenomena appearing in completely different 
contexts. We expect that this problem can be studied by finding the simplest model which 
describes the "essential nature" of the phenomena. Toward understanding the slugging phe- 
nomenon of granular flow, we analyzed phenomenological two fluid models in the previous 



papers [pl] , p2jl , and we derived the K-dV equation by focusing on the marginal point at which 
the uniform fluidized state becomes unstable. Unfortunately, however, phenomena which 
the K-dV equation can describe are quite restricted, and we were unable to explain slugging 
phenomenon even when taking into account higher order terms. Until now, a simple model 
describing slugging phenomenon has never been derived from phenomenological equations. 
Study of traffic congestion is in a similar situation. 

One of the reasons for such a poor understanding is that many of the phenomenological 
models to be analyzed are rather complicated. However, a quite simple model equation 
exhibiting traffic congestion has been recently proposed by Bando et.al. [|23|,24|j. In Section 
HI we introduce this model and review the study by Bando et.al. In Section III], we analyze 
this model equation and derive the modified K-dV (mKdV) equation at a critical point. This 
equation is universal in the sense that it is independent of the materials in question. However, 
since the mKdV equation is a completely integrable system, almost all solutions (solitons) are 
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structurally unstable, that is, they change discontinuously in response to perturbations of the 
system. Nevertheless, as shown in Section [IV], there are solutions which deform continuously 
when the mKdV equation is appended with higher order corrections. We call such solutions 
'deformable solitons', distinguishing them from the structurally unstable solutions. We 
show that there is a unique deformable kink soliton from which, by continuous deformation, 
a kink solution characterizing the essential features of congestion can be obtained. In the 
final section, we discuss several points, including the correspondence between our study and 
a slugging phenomenon in granular flow. 



II. MODEL 

In the present study, we employ a traffic flow model of the form 

Xi = a[U(x i+1 - Xi) - Xi], (1) 



which was proposed recently by Bando et al. [23,23]. Here, the dot represents differentiation 
with respect to time, i is the vehicle index, ordered according to vehicle position x^ (1 < i < 
N), and U is the velocity which the drivers prefer, depending on the inter- vehicle distance 
between the car in question and that immediately ahead of it, 6j(= x i+ i — Xj). Eq.([TJ) 
expresses that all drivers tend to adjust their vehicles' velocities to this velocity U with a 
constant sensitivity a. The function U is assumed to be monotonically increasing, with the 
supplemental conditions that U(oo) < oo and U(0) = 0. Although the theoretical analysis 
developed in this paper does not depend on the detailed form of U(b), we assume U(b) takes 
the form 

U(b) = tanhO - 2) + tanh(2), (2) 

when we need to discuss the system behavior concretely. Further we consider the situation 
that the vehicles run on a circuit of a length L, and focus on the "thermodynamic limit" 
that (L, N) — > (oo, oo) with L/N constant. Then, in the model system under consideration, 
there are two parameters, a sensitivity a and a mean distance b{= L/N). 



Eq.(|l|) has a uniform solution written as 

Xi(t) =bi + U(b)t + constant, (3) 

which expresses that all vehicles run with the same velocity U(b). Linear stability analysis 
of this uniform state was performed P3LE|]], and it was found that the solution (H) is linearly 



unstable when 

a < 2U'(b). (4) 



This unstable region in the parameter space (a, b) is shown in Fig.[3] |23|j24]| . Here, the prime 



refers to differentiation with respect to b. As seen in this graph, there exists a critical value 
of a, denoted by a c , such that uniform states with any b are linearly stable when a > a c , 
while uniform states for some b in a neighborhood of b c are unstable when a < a c . In 
the parameter space, the point (a c , b c ) is located at the peak of the marginal stability line 
defined by the equation a = 2U(b). We thus can determine the critical point by solving the 
equations U"(b c ) = and a c = 2U'(b c ). For later convenience, we present the expression of 
the relevant eigenvalue cr(k) of the m-th Fourier mode (k = 2^m/N) at the critical point 
(a c ,b c ) in the long wave-length limit (k — > 0): 

a(k) = °^ik + a -^e - - ^ + o(k% ( 5 ) 

System behavior under the condition that the uniform state is unstable has been inves- 
tigated through numerical simulations of Eq.(|l|) [23|,24}1. It has been observed in this case 



that vehicles separate into several regions in the circuit, in which vehicles form clusters with 



inter- vehicle distances fixed at the value bi or bu [p3| , p4| . Here, bi and bn satisfy &l < b < bu 
(see Figs. ^ and Since U{b L ) < U(bn), a cluster with an inter-vehicle distance b^ {bn) 
consists of slow (fast) moving vehicles. It has been interpreted that a cluster characterized 



by &£, suffers traffic congestion p3| , p4| . We note that each individual vehicle does not remain 
in a particular cluster. In fact, vehicles at the edge of clusters behave as follows. A vehicle at 
the head of a cluster with inter- vehicle distance bn encounters a slow-moving vehicle. Then, 



the driver decreases his vehicle velocity as the distance to the preceding vehicle decreases, 
and eventually this vehicle becomes a member of a cluster with inter- vehicle distance b^. 
Inversely, a vehicle at the head of a cluster characterized by b^ will become a member of a 
cluster characterized by bn after some period. That is, a given vehicle in the circuit experi- 
ences a series of slow- and fast-moving states. In this way, members of a cluster change in 
time. 

The most remarkable feature is that the inter- vehicle distances bi and bn are determined 



uniquely, irrespective of initial conditions, number of clusters, and b p3| , |24|| . In addition, 
there are two kinds of interfaces connecting two neighboring clusters, that is, one connecting 
the head of a bn cluster to the tail of a bi cluster, and the other connecting the head of a 
b L cluster to the tail of a b H cluster. It was concluded that 6^, bn and the shapes of the 
interfaces characterize traffic congestion p3[ , |24|| . In this paper, we regard the interfaces as 
kink and anti-kink solutions of Eq.(|l|) and develop a theory to determine &l, bn and the 
shapes of the kink. 



III. REDUCED EQUATIONS 

When we apply the methods of dynamical system theory to Eq. ([!]), it is convenient to 
rewrite it in the following form 

Vi = a[U(bi) -Vi), (6) 
bi = v i+1 - v h (7) 

where vi corresponds to the i— th vehicle velocity. We first note that the equality J2ibi = 
holds under periodic boundary conditions. Since this implies that J2i h is a conserved 
quantity, we expect that the characteristic time scale of the m-th Fourier component of the 
variable b goes to infinity as k —>■ ( see Eq. (|^) ) . On the contrary, from the form of Eq. ([]) , 
we conjecture that Vi relaxes to some value determined by the configurations of {bi} on a 
time scale of order I /a. Thus, when we are concerned with the long time behavior, we are 
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allowed to assume that the value of Vi is determined adiabatically by {hi}. This statement 
is expressed by the relation 



Substituting Eq.(|8|) into Eqs.(|J) and (|7|), we obtain the equations for V, and we can solve 
these equations by formal polynomial expansions in 6j. However, due to the complicated 
nature of these expressions, it is difficult to extract the essential nature of solutions. Thus 
we will focus on the system behavior near the critical point (a c , b c ) to simplify the problem. 
With such a simplification, the nature of kink solutions can be described, as we shall see 
later. Carrying out numerical simulations of Eq. ([!]), we found that a kink solution appears 
supercritically below the critical point a = a c . Therefore we can expect that each of the 
variables (bi,Vi,i,t) possesses a self-similar nature at the supercritical bifurcation point. We 
seek to determine the scaling properties of the variables when a = a c (l — e 2 ), where we have 
introduced e as a small scaling parameter. From the parabolic form of the marginal stability 
line around the critical point, it is expected that the amplitude of a kink solution is scaled 
as \b — b c \ oc e. Since the characteristic scales of space and time diverge at the critical point, 
using a(k) given by Eq.(|5]), we determine the scalings of i and t in the following way. First, 
the real part of Eq.([5]) contains both a fourth order dissipation term and a negative diffusion 
term whose coefficient is proportional to a — a c . From the balance of these two terms, k 
scales in proportion to e. This leads to the scaling relation i oc e . Second, the scaling 
of t is determined by the lowest order term in Eq.(||) , which is given by a dispersion term 
proportional to k 3 when we eliminate a propagation term by shifting to a moving coordinate 
system. Thus, t is scaled as t oc e~ 3 . 

Based on these considerations, we express the variable b{ as 



where U'{b c ) > 0, U"'{b c ) < 0, and Z and T are a scaled position and time defined by 



Vi{t) = V(bi, b i+ x, bi-t, ■ ■ ■)• 



(8) 




(9) 



Z = 2e(i - U\b c )t) 



(10) 



6 



T = -e 3 U'(b r )t. 



(11) 



Here, elaborate prefactors in Eqs.(PD,([TD|) and ( |TTD are introduced so that we can obtain the 
simplest form of the equation for b. (See Eq. (|l8D .) Then, from Eq.(U), Vi is expressed as 

Vi = V(e b(Z, T), e b(Z + 2e, T), e b(Z — 2e, T), • • ■). (12) 

Since e is a small parameter, a Taylor expansion can be applied to the terms like b[Z + 2e, T) 
in Eq.(|T2|). This leads to the expression 



Vi = V(eb(Z,T),e 2 d z b(Z,T),e 3 d 2 z b(Z,T) 



Substituting Eq . (|13"D into the right-handed side of Eq.(^), we obtain 



(13) 



e 4 d T b = A[V(eb(Z + 2e, T), e 2 d z b(Z + 2e, T), 
= ed z B(eb,e 2 d z b,e 3 d 2 z b,---), 



V(eb(Z,T),e 2 d z b(Z,T),---)], 



I, (14) 

where A is a constant given by A = (8U\b c ) / 3) U' (b c ) / \ U m (b c )\. Here, since V remains 
unknown, we determine by assuming the following form: 



V=Y / e n V^(b,d z b,d 2 z b,---). 

71=0 

From the form of Eq.(|T^), we find that the terms V^ n ' are expressed as 

= V , 

V« = Vj), 

1/(2) = V 20 b 2 + V 21 d z ~b, 

V® = v 30 b 3 + V 31 d z b 2 + V 32 d 2 ~b, 



(15) 



(16) 



where Vb, Vi, V20, V21, V30, V31, V32, etc. are constants which are calculated in the following 
way. First, we substitute Eq.flToD ( using Eq.([L6|) ) into Eq.fl^), where i)i is expressed as a 
functional of b, using Eqs.flTBD and (0). Then, collecting terms of equal order on both sides, 
noting the relation a = a c (l— e 2 ), and comparing coefficients of b, b 2 , d z b etc., successively, we 



obtain V^ n ' . ( Terms on both sides of the equation are now functionals of b. ) Substituting 
V derived in this way into Eq . (|14D , we get a reduced equation for b in the form: 



d T b = «9 Z £ e n B^(b, d z ~b, d 2 z ~b, ■■■). (17) 



n=0 



Since we know the scalings of (b, Z, T) , no explicit calculation is necessary to find which 
terms can appear in B^°\ Such terms must be of the same order as dxb, and therefore the 
three terms b ,dzb 2 and d\b can appear in B^°\ but the term dzb cannot. This can be 
explained in the following way. Since U(b) — U(b c ) is an odd function of b — b c at least up 
to the order e 3 at the critical point, V has the same symmetry up to the same order in e. 
This implies that V20 and V31, defined in Eq. ([16|) , vanish. As a result, also has the 
same symmetry at the critical point. An explicit calculation in fact does yield a B^ which 
contains only b and d\b . Substituting this B^ into Eq.(|T7]), we finally get the modified 
K-dV (mKdV) equation, 

d T ~b(Z, T) = d 3 z b(Z, T) - d z ~b 3 (Z, T), (18) 



which has been derived in the contexts of plasmas, or helical magnet systems, etc. |25|-|27 . 

We consider the steady traveling solutions of the mKdV equation (0). Substituting 
b(Z,T) = g(Z - cT) into ([18]), we get 



g-cg + 7, (19) 



dZ 2 

where 7 is an integral constant which is determined by the boundary conditions. Since the 
periodic boundary conditions we assume give 7 = 0, we consider the case 7 = in the 
argument below. Here, it is worth noting that the steady solutions of the non-conserved 
Ginzburg-Landau equation also satisfy Eq.(|T9D with 7 = 0. Also, Eq.(|T9|) can be regarded 
as the Hamilton equation for a nonlinear oscillator in a fourth order potential and can be 
solved analytically. When c < 0, the potential is convex everywhere, and thus all solutions 
are unbounded except that one which just remains on a peak. When c > 0, the potential 
has double peaks, and solutions bounded between these two peaks can exist. As a result, 
the following three parameter family of solutions ( with 7 = 0) are obtained. 
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g(Z - cT) = a «\e; m) = ,/-^ S n f ; m) , (20) 

V 1 + m \ Zir / 



where sn(x,m) is a Jacobi elliptic function, with sn(x, 0) = sm(x)andsn(x,l) = tanh(x), 
and K(m) is the complete elliptic integral of the first kind |28]. This family is parametrized 
by a traveling velocity c (c > 0), a modulus m (0 < m < 1), and a phase (0 < # < 27r). In 
particular, solutions with m = I take the form of kink patterns. Thus, there are an infinite 
number of kink solutions parameterized by the traveling velocity c and position </>. 

Now, let us recall that in numerical simulations of Eq. ([!]), the amplitude of kink solutions 
is determined uniquely irrespective of initial conditions. This leads us to the question: which 
solution correspond to the observed one? For such a selection problem, first of all, the 
stability of the solutions must be discussed. We note here that the mKdV equation is a 



completely integrable system ||29|| . That is, the solutions are "neutral", in the sense that 
they are specified by an infinite set of conserved quantities whose values are determined 
by the initial conditions. This implies that the stability argument for the solutions cannot 
identify the observed one. Then, it must be noted that integrable systems are structurally 
unstable. If higher order terms of the equation for b are taken into account, the phase space 
structure changes discontinuously. We are now led to consider the mKdV equation with 
higher order corrections. After some calculation, we obtain 

b = dfb- dzb 3 - -\d\(b + d 2 z b - ab 3 ) + (3d z b A } + 0(e 2 ), (22) 
a 

where a = 2/3 and (3 = when the an explicit form of U(b) given in Eq.(Q) is assumed. 
Here, a = 2/3 holds for an arbitrary form of U(b) satisfying U"'(b c ) < 0, but (3 is derived as 



(3 = (U""(b c )/3)\JU'(b c )/\U'"(b c )\ 3 . In the argument below, we consider the general case, in 
which P^Q. 
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IV. DEFORMABLE KINK SOLITON 



We now consider how the mKdV steady traveling solutions (|20|) are affected by the 
perturbation resulting in Eq. (f22|) . Substituting b(Z,T) = g(Z — cT) into Eq.(eqn:mkdv- 
dissipation) , we obtain the equation for the steady traveling solutions as 



d 2 g , e 

-T7^ = 9 -cg + - 
dZ 2 a 



9 + ^-*9 3 )+P9 4 



+ 7- (23) 



dZ y dZ 2 

As before, we will concentrate on the case 7 = 0. We solve Eq.(^) using a perturbation 
method which is similar to the Krylov-Bogoliubov-Mitropolsky method for nonlinear oscil- 



lations |30|. First, we consider this equation without the perturbation. Then, since the 
differential equation for g is second order, the solutions can be expressed as trajectories in 
a two dimensional phase space. For this purpose we must introduce a coordinate system in 
the phase space. While the obvious coordinate system is (g,dg/dZ), we are not limited to 
this one. Since we are free to choose coordinate systems, we choose that which is best suited 
to treat the present problem. In the coordinate system we use, g(Z) is expressed as 

g(Z)=gW(6(Zy,m(Z)), (24) 

where the function g^ is given by Eq. (|20|) . Using the coordinate system (9, m), we rewrite 
Eq.flTD as 

^ = ^°)(m), (25) 
dm . . 

Iz = a < 26 » 

Next, we take the perturbation into account. Since a third order derivative term is 
included in the perturbation, the dimension of the phase space becomes three. However, 
noting that the motion in the new dimension has a small time scale of order e, we can treat 
Eq.(p3|) as a two-dimensional system in our perturbation theory. That is, we express g as 

oo 

g = g(°\6;m) + Y,e n 9 in) (0;ml (27) 

n=l 
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and the evolution of the coordinates 9 and m are written as 

d6 



Q, = fi (0) (m) + Y, e ™ fi(n) ( m ) , ( 28 ) 

n=l 

dm 
dZ 



dZ " 

= M= ^e n M (n) (m). (29) 

n=l 

Here, fit") and M( n ) can be determined by Eq.(p3|) perturbatively. We now present the 
method to calculate g w , and M {1 \ 

Substituting Eq.(|27|) into Eq.(|23|) and collecting terms of order e, we obtain 

LgV>=hP\ (30) 

where 

L = -nW 2 d 2 e + 3gM 2 -c, (31) 

and 

= |2fi(°) (0«a| + M^d e d m ) + (d m n<®) M^d e - i (l + fiW 2 ^ - 3a</<°> 2 ) 9, + ^ 

(32) 

Here, from the fact that (9,m) parameterizes the family of solutions of Eq. flT9"|) , we can 
easily derive the relations L ■ dgg^ = and L ■ d m g^ = 0. This implies that there is a 
solution only under a solvability condition. Since we are interested in the case that the 
perturbation expansion is well-defined, we enforce the solvability condition expressed as 

{a B gW,hW) = {a n gM,hW) = 0, (33) 

where (/, g) = fjf* d9f(9)g(9) for 2ir periodic functions / and g. These two equations 
yield the expressions of and M«. Then, can be obtained from Eq. (|30"D , owing to 
the solvability condition. We show only the result of MW which will be important in the 
following discussion. We obtain 

= r(m)[c*(m) — c], (34) 
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where 



m 



r(m) 



(1 + m) 5 /2 J 



IK 

dx 



(cPsn(x; m)) 2 + 6ma(d x sn(x; m)) 2 (sn(x; m)Y 



ad„ 



m 



(1 +m) 3 / 2 Jo 



■IK 



dx(d x sn(x; m)Y 



(35) 



and 



c m 



1 + m) / cfxfc^snfx; m)V 

jo 



4X 



(<9 2 sn(a;; m)) 2 + 6ma(<9 x sn(x; m)) 2 (sn(x; m)) z 



(36) 



The profiles of the functions c*(m) and r(m) are displayed in Fig.[|. For later convenience, 
we define m*(c) as the inverse function of c*(m). 

We study the differential equation Eq.(p9|) with Eq.(^) under the boundary condition 
miZ = 0) = mo- We discuss the behavior of solutions for four cases: (i) there is no m*(c) 
for a given c, (ii) m = m*(c), (iii) m < m*(c) and (iv) m > m*(c). For case (i), the 
sign of M^(m(Z)) is positive (negative) definite for all Z. When M^(mo) is positive 
(or negative), m(Z) increases as Z increases (decreases). Then m(Z) = 1 is realized at a 
finite Z. Since (m(Z)) is positive (negative) even at this point, m(Z) overshoots 1 as 
Z increases (decreases) further. Thus a perturbed solution cannot be defined, and there 
is therefore no steady traveling solution in this case. For case (ii), M'^mo) = 0. This 
leads to the equality m(Z) = m*(c) for all Z, that is, the perturbation does not make m 
inhomogeneous. Thus, the expression 



,(°) 



(9,m*(c)) + eg (1 \9,m*(c)), 



(37) 



with 9 = Q(Z — cT) + gives a steady traveling solution of Eq. (|2"2"|) , where 



Q = fi(°)(m*(c)) +eQ w {m*{c)). 



(38) 



We call the unperturbed solution with (c, m*(c)) a "deformable soliton" because the effects 
of infinitely small perturbations can be absorbed into infinitely small deformations of the 
unperturbed solution. The corresponding perturbed solution is called a "deformed soliton" . 
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In case (iii), M^(m ) is negative. Here, m(Z) increases (decreases) as Z decreases (in- 
creases) and converges asymptotically to the fixed point m(Z) = m*(c) (m(Z) = 0). Thus, 
we obtain m(—oo) = m*(c) and m(oo) = 0, that is, the corresponding solution is a propagat- 
ing front solution connecting the solution characterized by (c, m*(c)) and the uniform state 
characterized by m = 0. In case (iv), (mo) is positive. In this case, m(Z) decreases 
as Z decreases and converges asymptotically to the fixed point m(Z) = m*(c) similarly to 
case (iii). Nevertheless, m(Z) increases as Z increases and overshoots 1 similarly to case (i). 
Thus a perturbed solution cannot be defined, and thus there is no steady traveling solution 
in this case, as in case (i). As a result, all solitons except deformable solitons are destroyed 
by the infinitesimal perturbation. The deformable solitons of the mKdV equation form a 
two parameter family parametrized by the modulus m and arbitrary phase 0. The traveling 
velocity of the soliton with modulus m is given by c*(m). Note that the soliton with m = 1 
represents a kink pattern, while the other solitons are periodic. Therefore, we can expect 
that the deformed kink soliton corresponds to that observed in numerical simulations of the 
traffic flow model, Eq.fll]). 

In order to support this conjecture further, we need to check the linear stability of the 



deformed solitons by studying Eq.(p2). However, we have not yet obtained a rigorous result 
for this problem. In this paragraph we present a plausible conjecture. First, one may 
prove the asymptotic stability of all kink solutions of the mKdV equation by introducing 
a physically acceptable norm, as recently done for solitary pulse solutions ]3"I| . Then, the 
deformed kink soliton, which is a solution of Eq. (|22|) , can be shown to be stable if it can be 
proven that the correction of eigenvalues associated with the linear stability is bounded to 
order e. On the other hand, from the fact that the tail of the kink solution simply decays 
exponentially, we expect that the kink interacts with a neighboring anti-kink attractively 
||32|| . This implys that periodic deformed solitons with sufficiently large m are unstable. 
Then, since there is no apparent reason to assume that there is a critical m below which the 
stability of these solutions is changed, we conjecture that all deformed periodic solitons are 
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unstable, while the deformed kink soliton is stable. 

Next we confirm that the deformed kink soliton corresponds to that observed in numerical 
simulations of Eq.(|l|). First, using Eq. (p0|) , we find that the profile of hi corresponding to the 
deformable kink soliton is a kink pattern connecting the two uniform states characterized 
by and where bi = b c — Ab, and b H = b c + Ab with 



Ab = 2e, 



c*(l)U'(b c ) 

~~\\ \u'"{l c )\ ■ {6J) 

This expression can be compared with that obtained by solving the model Eq.(|I|). We thus 
carried out numerical simulations of Eq.(|l|) using the fourth order Runge-Kutta method. In 
Figs.[| and || the result is shown together with the theoretical curve. It is seen that Eq. (p9|) 
gives a good approximation of the numerical result for e smaller than about 0.25. Although 
the discrepancy increases as e becomes large, we believe that the behavior for larger e can be 
accounted for simply by calculating higher order terms in Eq . (|T7|) . Therefore, we conclude 
that the kink pattern arising in traffic congestion corresponds to the deformed kink soliton. 

V. DISCUSSION 

To this point, we have analyzed Eq.(^) around the critical point (a c ,b c ). Similarly, we 
can obtain the result that system behavior near another point on the marginal line defined 
by Eq.([y]) is described to lowest order in a perturbative expansion by the K-dV equation 
|pl| , |22|j33| . The K-dV equation is also a completely integrable system with an infinite number 
of conserved quantities f34] , |35|| . Then, owing to the effects introduced by higher order terms 



in the expansion, a pulse solution with some amplitude is thought to determine the long 



time behavior |pq , p7|l . Here, we should remark that this K-dV pulse solution is observed 
only in a quite narrow region in parameter space. This contrasts sharply with the case of 
our deformed kink soliton, which is observed in a wider region, even beyond the marginal 
stability line. ( This is discussed further below. ) We believe that understanding the 
difference between these cases is important. However, we have no theory related to this 
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problem. We conjecture that a crucial step is to formulate the non-perturbative nature of 
continuous deformations of the solution. 

Further, as shown in Figs.[7| and |8|, kink patterns are observed even when the uniform 
state is stable. This result implies that final states depend on initial conditions, because 
the non-uniformity of vehicle velocities vanishes when the magnitude of fluctuations around 
the uniform states are sufficiently small. Moreover, if a kink pattern forms, its amplitude is 
determined by only the value of a irrespective of the linear stability of the uniform states. 
Then, from Eq.fl39f), we can expect that the parameter region where kink patterns can appear 
is approximated well as 



recalling that a = a c (l — e 2 ). The validity of this expression was confirmed by numerical 
simulations. On the other hand, the marginal stability line of uniform states is given by 



kink patterns can appear even when the uniform states are stable. 

Finally, we briefly discuss a slugging phenomenon appearing in quasi-one dimensional 
granular flow. As a typical experimental system, we consider the sedimentation of solid 
particles (powder) in fluid confined in a narrow pipe. We assume that the density and 
velocity of solid particles are dynamical variables which can be regarded as a functions of 
the vertical position z owing to the system's quasi-one dimensional geometry. In order to 
see the correspondence with the traffic flow model further, it is convenient to introduce a 
description employing "powder elements" instead of a continuous description. That is, the 
velocity of powder v(z, t) contained in the i-th element located at the position z is represented 
by Vi(t). (This is nothing but a shift from the Euler picture to the Lagrange picture.) Let 
us consider the equation of motion for the powder element. The most important force, 
responsible for the variety of dynamical behavior, is the drag force due to the fluid around 




(40) 



b-b, 




(41) 



This line is contained in the region defined by Eq.(^Dj) because 2c*(l) = 5/2. Therefore, 
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the settling powder element. This is proportional to the velocity difference between the 
powder element and fluid. Further, it is known that the proportionality coefficient (drag 
coefficient) is determined by the powder density. Taking the drag force and the gravitational 
force into account, we obtain a simple form for the equation of motion: 

v% = -a s {p)Vi - g = a s {p)[U s {p) - (42) 

where a s (p) is a (rescaled) drag coefficient, g is a specific gravitational acceleration, and 
U s (p) = —g/a s (p). This model can describe a uniform state in which all elements settle 
with the same velocity U s (p). There is one unknown function a s (p) in the model, but this 
can be determined experimentally. Actually, from experiments on sedimentation of solid 
particles in three dimensional fluid f39| , it is known that the sedimentation velocity U s is a 



monotonically decreasing function of the powder density. Using this experimental data, the 
drag coefficient a s (p) can be determined. 

One may suspect that the model equation ( |4"^) is too simple to describe complicated 
phenomena and that for a realistic description we must include a dissipation force accounting 
for collisions between particles and a chemical potential force. Still, our model equation ( f4"2"D 
is instructive when we consider the correspondence between slugging phenomenon in granular 
flow and traffic congestion, since Eq.(^) is of the same form as the traffic flow model (|l|) if 
we neglect the p dependence of a s . (Note that p may be related to inter- vehicle distance b.) 
The problem which needs to be considered is whether or not the analysis developed above 
can be applied to the model system (0). The main difference between the two models 
is the p dependence of a s , which causes the asymmetric term dzb in B^°\ ( Recall the 
argument appearing before Eq.([L8j). ) Due to the presence of this term, we cannot obtain 
the mKdV equation even when we focus on a critical point. Nevertheless, we conjecture 
that the asymmetric term does not change the essential nature of solutions, i.e., the effects 
caused by the asymmetric term can be absorbed into deformation of the kink soliton in the 
mKdV equation. For this reason, we believe that slug flow is described by a deformed kink 
soliton. In order to prove this statement, we need to develop a non-perturbative method 
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describing the continuous deformation of solutions. To construct such a theory is a future 
problem. 
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FIGURES 

FIG. 1. Marginal stability line of uniform states in parameter space. The critical point (a c , b c ) 
is marked by the filled box. 

FIG. 2. Spatio-temporal pattern exhibited by Eq.(|]), with a = 1, L = 70, and N = 35. 
Positions of vehicles at every 5 time steps are indicated by dots. The parameters are chosen such 
that a uniform state is linearly unstable. 

FIG. 3. Profile of vehicle configuration when t = 1000 (circles with broken line) with the same 
parameters as in Figj2|. For the sake of reference, initial conditions (diamonds with solid line) are 
also displayed. The horizontal axis i and vertical axis b\ represent the index of vehicles and the 
distance between the i-th and (i + l)-th vehicle, respectively. 

FIG. 4. c*(m) in Eq.(|36|) versus m (solid line) and r(m) in Eq.(|35]) versus m (dashed line). 
Vertical axes of solid and dashed lines are graduated on the left and right, respectively. Note that 
c*(l) = 5/4. 

FIG. 5. Amplitudes (max|& — b c \) of the kink patterns observed in numerical simulations are 
plotted for e (diamond). The solid line shows the theoretical curve given by Eq.(|39|), 

FIG. 6. Profile of kink pattern observed in numerical simulation with e = 2~ 3 , N = 128, and 
L = 256 (circle). The solid line shows the deformable kink soliton obtained theoretically, whose 
phase (i) is adjusted to match that of the numerical one. 

FIG. 7. Spatio-temporal pattern exhibited by Eq.([l|), with a = 1,L = 70, and iV = 64. 
Positions of vehicles at every 5 time steps are indicated by dots. The parameters are chosen such 
that a uniform state is linearly stable. 
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FIG. 8. Profile of vehicle configuration when t = 1000 (circles with broken line) with the same 
parameters as in Fig]7|. For sake of the reference, initial conditions (diamonds with solid line) are 
also displayed. The horizontal axis i and vertical axis b{ represent the index of vehicles and the 
distance between the i-th and (i + l)-th vehicle, respectively. 
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T.S.Komatsu and S.Sasa Figure 1. 
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